Excitonic optical spectrum of semiconductors obtained by time-dependent 
density-functional theory with the exact-exchange kernel 
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Applying the novel exact-exchange (EXX) Kohn-Sham method within time-dependent density- 
functional theory, we obtained the optical absorption spectrum of bulk silicon in good agreement 
with experiments including excitonic features. Analysis of the EXX kernel shows that inclusion of 
the Coulomb coupling of electron-hole pairs and the correct long-wavelength behavior in the kernel 
is crucial for the proper description of excitonic effects in semiconductors. 
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The calculation of optical spectra of solids has been 
traditionally one of the most important and challenging 
areas of first-principles material investigations. The in- 
teraction of excited electrons and holes plays a crucial 
role in optical excitations, however properly incorporat- 
ing such effects ab initio is theoretically and computa- 
tionally very difficult. In the past several years there has 
been a major advance in the field W based on the so- 
lution of the Bethe-Salpeter equation (BSE) 0] for the 
two-particle Green's function starting from quasiparticles 
obtained within the GW approximation (GWA) ||, but 
unfortunately the calculation scheme is by nature com- 
plicated and demanding. 

One attractive alternative is extending density- 
functional theory (DFT) to the time-dependent (TD) 
case j|. TDDFT is a well-founded theory for the treat- 
ment of electronic excitations in general and, unlike 
the GWA-BSE route which is presently implemented as 
a post-DFT method, provides a computational scheme 
based directly on Kohn-Sham (KS) one-particle states. 
Moreover, once the KS one-particle states are given, 
the electronic linear response properties are completely 
determined by the orbital-independent Hartree kernel 
Fn(r, r') = l/|r — r'| || and the dynamic exchange- 
correlation kernel 



F xc {r, r'-t-t') 



Sv xc (r;t) 
Sn{r';t') ' 



(1) 



Thus, in principle, one only needs to find a good ap- 
proximation to the exchange-correlation kernel for prac- 
tical TDDFT applications, which represents a significant 
conceptual and computational simplification. In prac- 
tice, however, the conventional adiabatic local-density 
approximation (LDA) and generalized-gradient approx- 
imation (GGA) kernels are known to be inadequate for 
the study of optical excitations in solids. A primary ex- 
ample is their incorrect long-wavelength behavior for in- 
sulators which has been recently emphasized in several 
contexts g g, §. 

Deficiencies of the LDA and GGA in fact appear 
already on the level of the ground-state exchange- 
correlation energy E xc and potential v xc (r) = 



SE XC /Sn(r). With respect to energetics, for example, the 
LDA and GGA inherently fail to describe the quasi-two- 
dimensional electron gas || and the van der Waals inter- 
actions [[l0| . For v xc and the corresponding KS eigenvalue 
spectrum, LDA and GGA v xc incorrectly decay exponen- 
tially rather than as — 1/r in finite systems and accord- 
ingly their highest occupied orbital energies are too high 
and unoccupied orbitals do not exhibit Rydberg series, 
and for solids their band gaps are too small. The de- 
ficiencies of the LDA and GGA v xc in particular pose 
significant practical difficulties for the study of electronic 
excitations within the DFT formalism. 

In this regard, recent realizations of multi-dimensional 
KS exact-exchange (EXX) 0, |f, || and effective 
EXX 0, [HI [u| [l7| methods provide an interest- 
ing opportunity. The orbital-based self-interaction-free 
EXX methods not only provide realistic local multiplica- 
tive KS exchange potentials and KS eigenvalue spectra 
for molecules @, [if, H but also 

give band 

gaps of semiconductors in good agreement with experi- 
ments |ll], |l4| . We have also recently shown [pi that the 
EXX spectrum at the one-particle level without any pre- 
viously applied post-DFT modification such as a quasi- 
particle shift pl| gives a very good description of the ab- 
sorption spectrum of semiconductors excluding excitonic 
features resulting from two-particle interactions, which 
supports the notion that KS eigenvalue differences rep- 
resent well-defined approximations for excitation ener- 
gies pfj| , However, the remaining excitonic character 
in the spectrum is not properly treated by the LDA or 
GGA kernel, so a complete set of DFT methods for the 
investigation of electronic optical excitations in solids is 
still missing. 

Given the importance of proper inclusion of excitonic 
effects in the study of electronic optical excitations and 
the encouraging performance of the EXX method for 
band structures of semiconductors, we consider in this 
work a TDDFT scheme based on the EXX kernel. At 
the formal level, in contrast to LDA and GGA kernels, 
the EXX kernel is nonlocal in real space, depends explic- 
itly on the frequency, and has the correct long-wavelength 
limit behavior |?], |22|. Taking bulk silicon as the repre- 
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sentative semiconductor system, we show that the EXX 
kernel indeed provides optical spectra of semiconductors 
in good agreement with experiments including excitonic 
effects and further analyze the origin of the excitonic peak 
and the validity of locality approximations in space and 
time. We additionally make comparisons between EXX- 
TDEXX and GWA-BSE, which renders new insight into 
the TDDFT approach in general Q. 

To obtain the optical spectrum of a solid, we computed 
the frequency w-dependent macroscopic dielectric func- 
tion £m(w) which can be written in terms of the modified 
full linear response matrix £ as 



1 



lim F H (q)x(q;^)|G=G' 

9^0 



(2) 



where q is the photon momentum, and -FH(q) = 
Fh(G, G',q) = <5 G , G '47r/|q+ G| 2 . Within TDDFT, % 
is completely determined by the KS linear response ma- 
trix Xo and the kernel matrix F xc according to 



X(q;w)= 1 - Xo(q;w){^H + F xc (q;cj)} xo(q;w), 



(3) 



where F H = for G = G' = and F H = F H oth- 
erwise. At the independent-particle level, \ = XO) the 
two-particle interaction effects due to -Fh and F xc are 
both ignored, while only the F xc part is ignored at the 
time-dependent Hartree level. As is apparent from Eqs. 
(H) and (j^) , the long- wavelength q — > behavior of the 
"head" (G = G' = 0) and "wing" (G = and G' ^ or 
vice versa) elements of F KC (in insulators, 0(l/q 2 ) for the 
head and 0(l/q) for the wings [0, ^2| ) plays an important 
role. 



The computation of xo is straightforwardly done with 
KS one-particle states and energies, and consequently 
F xc remains as the only important quantity to be deter- 
mined. So far, for F xc , the adiabatic LDA kernel F XC DA = 
S 2 F^ A I (Sn Sn') has been almost exclusively employed. 
It is local in real space and frequency- independent, which 
results in a reciprocal-space representation independent 
of q and uj, F^ C DA {G, G', q; u) = F£ C DA (G - G'). While 
numerical advantageous, this drastic simplification is 
known to be problematic for the study of optical re- 
sponses in solids. For example, the head and wings of 
the LDA kernel are incorrectly nondivergent in the q — > 
limit and this defect cannot be corrected by a semilocal 
GGA kernel. 

Instead, we adopted in this work the EXX kernel, 



F™(G,G',q; W )= ]T ^(G.Gi.qjo;) 

G!,G 2 (4) 

x fl- x (Gi, G 2 , q; w) Xo\G 2 , G', q; w), 

which is nonlocal in real space and explicitly depends 
on the frequency. The expression of the EXX kernel 
"core" 7J X , which was interpreted in the many-body 
diagrammatic picture as the first-order self-energy and 
vertex corrections to the irreducible polarizability x — 
Xo + Xo F xc x Q j has been presented in Ref. [^2j for the 
case of periodic insulators. We have further shown that 
F x xx has the q — > behavior of the exact F xc and thus 
rectifies one serious deficiency of the LDA and GGA ker- 
nels. The full expression for H x is rather lengthy and we 
list here only the three resonant terms : 



^-"-(qjw) 
H x B - res (q; W ) 



EE[ 



(ak|e- 4 ( q+G )- r |sk + q) (sk + q; bkf\w c \tk' + q; ok) (tk! + q|e i ( c '+ G ')- r |6k / ) 



sk btk' 



(e a k - £sk+ q + + i5)(ebk< - e*k' 



i5) 



(ak|e- 4 ( q+G )- r |sk + q)(bk\t x - v x \ak)(sk + q|e l ( q+G '> r |6k) 



afcsk 



(e Q k - e s k+ q + uj + i8)(ebk - e s k+q - 



iS) 



(5) 



n 

asik 



(ak|e- 4 ( q+G )- r |sk + q)(sk + q\E x - %\tk + q)(tk + q|e l ( q+G ')- r |ak) 



(e Q k - e s k+ q + w + iS){e a k - e*k+ q + w + id) 



where two is the spin factor, Q is the crystal volume, 
{a, b} are valence bands, {s, t} are conduction bands, the 
matrix elements (i; j\wc\l] m) are four-index integrals de- 



fined as 

(ik + q;jk'\w c \lk' + q; mk) = Jdrf dr' 0* k+q (r)^* k ,(r') 

x wc(r, r')</>ik'+ q (r)0mk(r')> 
(6) 

Yi x denotes an orbital-dependent exchange operator of 
the form of the Hartree- Fock-like exchange operator but 
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constructed from the KS orbitals, i.e., 



(ik + qlS^ljk + q) = - ^(ik + q; ak'\w c \ak': jk + q) 



ok' 



(7) 



and v x is generated by the orbital-independent EXX KS 
potential v x (r). Here «ic(r,r') yielding wc denotes the 
generalized Coulomb interaction, e.g., l/|r — r'| for the 
bare Coulomb interaction. 

Compared with the standard BSE approach, although 
we also have to perform the four-index Coulomb inte- 
grals of Eqs. (§) and (0), we can calculate optical spectra 
without diagonalizations in the space of occupied and un- 
occupied single-particle states. Due to the huge number 
of k-points involved in the optical spectrum calculations, 
the size of the matrix to be diagonalized can be very 
large, and the fact that we can avoid the diagonalization 
process altogether in principle represents a significant nu- 
merical advantage. In addition, note that F^ xx is free 
of Coulomb singularities |22j and thus the formal validity 
of the EXX TDDFT approach does not depend on the 
consideration of the thermodynamic limit. In the BSE 
equation, on the other hand, integrable singularities ap- 
pear which require special numerical care (l| . 

We have implemented the full dynamic EXX kernel in 
reciprocal space employing the plane wave basis. The 
accuracy of the code has been carefully checked by nu- 
merically testing if the calculated exchange kernel acts 
as a functional derivative of the exchange potential with 
respect to the electron density at uj = 0, 8v^ xx = 
F^ xx Sn. To obtain F^ xx wc first carried out self- 
consistent EXX ground-state calculations at the experi- 
mental lattice constant of Si, 5.43A, and generated the 
KS potential p, 18 1. Ten special-k points and an orbital 
kinetic energy cutoff of 12.5 Ha have been used. In the re- 
sponse calculation step, we solved the KS equations once 
more at a larger number of k-points and obtained KS 
orbitals and eigenvalues as the input for the construction 
of xo and F^ xx . We adopted a shifted uniform k-mesh, 
for which we employed up to 22 x 22 x 22 k-points for 
Xo and up to 9 x 9 x 9 k-points for F^ xx . The kinetic 
energy cutoff of 10 Ha and 8 = 0.15eV have been used, 
and ten conduction bands were included. 

The first question we addressed with small k-points 
sets (e.g. 5x5x5) was whether the adiabatic approxi- 
mation, defined as the u> = limit of Eqs. (||) and (||), 
is justified. Although the nature of the approximation 
within TDEXX is different from that of the BSE ap- 
proach in which the static screened interaction is taken, 
we found that the adiabatic EXX kernel generates an 
overall similar spectrum as the dynamic one. Computing 
F® only at a single u> results in a significant reduction 
of computational workload, and thus the converged cal- 
culation was performed in the adiabatic approximation. 

Figure [I] shows the absorption spectrum of Si ob- 
tained at the single-particle EXX level and by tak- 



ing into account two-particle interaction effects via Fr- 
plus F^ C DA (EXX+TDLDA) and via F H plus F^ xx 
(EXX+TDEXX). Compared with the much discussed 
LDA (and LDA+TDLDA) spectrum, which is incorrectly 
shifted to the lower frequency region by about 1 eV 
due to its well-known band gap underestimation, the 
single-particle EXX absorption edge and second (E 2 ) 
peak are in excellent agreement with the experimental 
curve Q due to the realistic EXX band structure [ pj| . 
However, the first (E\) peak originating from electron- 
hole attractions is much underestimated at the single- 
particle EXX level, and it is not recovered by taking into 
account electron-hole interactions via TDLDA, demon- 
strating the failure of F XC DA to describe excitonic effects. 
The EXX+TDEXX spectrum, on the other hand, is in 
excellent agreement with the measured data: The ab- 
sorption strength at the E\ peak region is correctly en- 
hanced while that of the higher frequency region is re- 
duced, which shows that F® xx indeed provides the cor- 
rect description of the important excitonic effects. 

To obtain the EXX+TDEXX spectrum in Fig. |, we 
used a slightly modified bare Coulomb interaction as wc- 
Employing the bare Coulomb interaction resulted in a 
collapse of the spectrum due to too strong long-range 
Coulomb interaction of electron-hole pairs at different 
k-points. We therefore cut off the long-range Coulomb 
interaction for these pairs, or set to zero the contributions 
to the integrals of Eq. ( J6|) fo r which (G+k— k') lies within 
the first Brillouin zone pTj. Devising a more systematic 
strategy might be desirable in the future. 

Now, to understand the underlying mechanism of the 
encouraging EXX+TDEXX result, we analyze F^ xx by 
investigating the role of different contributions to H x as 
shown in Fig. |^. First, by taking only , we obtain 
a spectrum resembling that of the one-particle EXX but 
shifted upward by about 0.3 eV. No excitonic feature 
arises at this level. This spectrum is similar to that 
of the GWA, which is in agreement with previous pre- 
dictions [[|| ID that GWA results should be recovered 
with the ff,f part of the screened EXX kernel. How- 
ever, note that we have used an only slightly modified 
bare Coulomb interaction, and, while the GWA spec- 
trum results from independent quasiparticles that in- 
volve (iV± l)-electron excitations, the DFT "quasiparti- 
cles" spectrum already represents electron-number con- 
serving (A^-electron) one-particle excitations at the time- 
independent level |(], pi]] . Second, if we take only 
H A , the spectrum is strongly enhanced at the lower fre- 
quency region (< 4 eV) while it is much reduced at the 
higher frequency region (> 4 eV). The excitonic E\ peak 
appears predominantly while the E 2 peak is slightly red- 
shifted by about 0.2 eV. This shows that the Coulomb 
coupling of electron-hole pairs is precisely the origin of 
the excitonic E\ peak. Finally, we set the head and wings 
of F® xx to zero and thus make them q-independent as 
in the case of the LDA or GGA. By doing so, we obtain 
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a spectrum similar to that of the EXX+TDLDA without 
excitonic peaks or whatsoever. This demonstrates that 
taking into account the correct q — > behavior of F xc is 
crucial. 

We finally comment on other recent TDDFT works 
that have also obtained optical spectra of semiconductors 
in good agreement with experiments. First, de Boeij et 
al. have left the domain of conventional DFT and re- 
sorted to TD current-DFT [§§. Albeit they had to in- 
troduce an empirical energy shift and a prefactor for the 
polarization functional, their approach may represent a 
different route to include the space nonlocality discussed 
above via a macroscopic functional. However, our result 
shows that employing the current response derived from 
a macroscopic functional is not necessary if one prop- 
erly incorporates the space nonlocality in the microscopic 
density response. In this context, recent work of Rein- 
ing et al. which was based on the idea of mapping the 
BSE to the TDDFT is in agreement with ours gj. But 
in their scheme they first invoked the GWA quasiparti- 
cle shift and second employed an empirical static scaled 
Coulomb kernel as F xc . Although their approach is not 
rigorous, we note that their simple static kernel never- 
theless exhibits a correct q — > divergence. 

In summary, we have reported calculations of the op- 
tical absorption spectrum of Si with the novel TDEXX 
approach. Of particular interest was the nature of the 
dynamic and nonlocal orbital-based EXX kernel which 
exhibits the q — > behavior of the exact F xc and success- 
fully generated excitonic features in the spectrum. We 
showed that including Coulomb coupling of KS electron- 
hole pairs at different k-points and the correct long- 
wavelength behavior of the exchange-correlation kernel 
is crucial, while the adiabatic approximation can be jus- 
tified. 
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FIG. 1: Calculated optical absorption spectrum of Si from EXX (dashed lines), EXX+TDLDA (dot-dashed lines), and 
EXX+TDEXX (solid lines) . Circles denote experimental data of Ref. |26[ 
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FIG. 2: Calculated optical absorption spectrum of Si from EXX+TDEXX with full EXX kernel (solid lines), only (dashed 
lines), only (dot-dashed lines), and full EXX kernel with head and wings set to zero (long-dashed lines). Circles: experiment. 



